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A  FULLY  TWO-DIMENSIONAL  EQUILIBRIUM  AND  TRANSPORT  MODEL  OF  THE 

POLO ID AL  DIVERTOR 


I .  INTRODUCTION 

Tokamak  designs  are  becoming  very  complex.  In  order  to  improve 

the  confinement  time  and  increase  the  beta,  present  and  future- 

generation  toroidal  confinement  devices  may  have  non-circular  cross- 

sections  or  multiple  magnetic  axes  or  both.  Vertically  elongated 

cross-sections  promise  improved  stability  coupled  with  higher  values 

of  the  plasma  current  and  total  beta.'*’  There  are  indications  that 

substantially  higher  plasma  pressures  can  be  more  confined  in  doublet 

2 

geometries  than  m  circular  geometries.  To  control  the  level  of 
impurities  and  thus  increase  confinement  time,  a  divertor  will  be 
needed  in  tokamak  reactors.3  These  essential  features  of  research 
devices  demonstrate  the  need  for  a  multidimensional  model  to  investi¬ 
gate  equilibrium  and  transport  in  present  and  proposed  experiments. 

The  simple  low-beta,  one-dimensional,  circular  flux  surface  model  is 
not  adequate  to  describe  the  physics  of  these  sophisticated  devices. 

In  an  elongated  or  "D"  shaped  discharge,  the  curvature  varies 
greatly  and  the  transport  coefficients  cannot  be  approximated  by 
surface  averaged  quantities,  as  in  the  usual  one -dimensional  schemes. 
In  a  high  beta  discharge,  the  flux  surfaces  change  shape  as  a  function 
of  time  and  a  multi-dimensional  scheme  is  required  to  model  this 
Manuscript  submitted  February  17, 1981. 
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effect.  The  doublet  and  poloidal  divertor  systems  contain  multiple 
magnetic  axes  with  a  separatrix.  They  are  inherently  two-dimensional. 


We  present  here  the  methods  we  have  developed  to  calculate  equili¬ 
bria  and  transport  in  a  fully  two-dimensional  geometry.  The  model  is 
sufficiently  general  that  it  is  capable  of  calculating  equilibria  and 
transport  given  rather  general  plasma  geometries  and  coil  structures. 

Here  we  apply  the  model  to  the  poloidal  divertor  system  (the  Princeton 
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Poloidal  Divertor  Experiment  (PDX)  )  and  describe  the  results. 

We  have  developed  a  fully  two-dimensional  Eulerian-Lagragian 
computer  simulation  model  of  tokamak  discharges.  This  computer  model 
is  used  to  investigate  the  time  evolution  of  a  finite-conductivity 
plasma  in  a  finite-conductivity  vessel.  It  uses  as  a  finite- difference 
mesh  a  general-connectivity  triangular  grid.  This  permits  the  inves¬ 
tigation  of  problems  of  much  greater  complexity  than  more  conventional 

one-  or  two-dimensional  schemes.  This  code  is  an  extension  of  the 

5  6 

previously  reported  Circular  tokamak  models  '  which  were  an  out- 

7 

growth  of  the  Linus  projects. 

The  code  we  have  developed  is  designed  to  follow  simultaneously  the 
resistive  diffusion,  transport  and  gross  dynamics  of  a  magnetically  con¬ 
fined  plasma  of  arbitrary  shape.  The  finite  resistivity  of  the  confining 
shell  permits  us  to  examine  the  discharge  for  times  exceeding  the  pene¬ 
tration  time  of  the  magnetic  fields  through  the  copper  shell.  The  model 
is  based  on  the  quasi-static  evolution  of  force  balance  in  the  plasma. 

The  simplest  set  of  equations  which  describe  a  plasma  in 
equilibrium  with  a  magnetic  field  is 
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(1-1) 


=  c  i x  I* 

_  4n  , 

V  X  B  =  —  2'  (1-2) 

V  •  B  =  0.  (1-3) 

Eq.  1-1  stems  from  the  equation  of  motion.  The  force  balance  between 
the  plasma  pressure  and  magnetic  forces  accurately  represents  a  system 
that  evolves  in  time,  provided  the  evolution  is  slow  and  nonideal 
forces  are  weak.  Generally,  this  implies  that  the  system  changes  on  a 
time  scale  characterized  by  the  resistive  diffusion  time.  The  question 
of  time  scales  is  treated  in  more  detail  in  Section  IV. 

The  usual  method  for  determining  this  equilibrium  consists  of 
casting  Eqs.  (1-1)  -  (1-2)  into  the  form  of  a  nonlinear,  elliptic,  par¬ 
tial  differential  equation  for  the  flux  function  ip.8  This  resulting 
equation  is  then  differenced  on  an  Eulerian  mesh  and  solved  by  itera¬ 
tion  subject  to  appropriate  boundary  conditions. 

Our  method  for  calculating  the  equilibrium  is  quite  different. 

The  equilibrium  equations  describe  the  relations  between  nested  toroidal 
surfaces9  (flux  surfaces) .  The  plasma  pressure  is  constant  on  these 
surfaces  and  B  and  j  lie  on  these  surfaces.  The  solution  to  the  equili¬ 
brium  problem  is  found  by  constructing  a  finite  set  of  flux  surfaces, 
each  denoted  by  a  flux  function  and  adjusting  those  surfaces  until 
they  are  in  equilibrium. 
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These  flux  surfaces  are  represented  by  line  segments,  whose  junc¬ 


tions  are  interconnected  to  construct  a  triangular  grid.  The  triangular 
gridding  scheme  allows  the  coordinate  system  to  be  tailored  to  the  prob¬ 
lem  being  investigated.  Complex  geometries  and  interfaces,  multiple 
magnetic  axes,  separatrices  and  limiters  (conductors  intersecting  flux 
surfaces)  are  readily  accommodated  with  this  model.  It  is  a  trivial 
matter  to  vary  the  resolution  across  tne  grid  and  thus  provide  high 
resolution  only  where  it  is  needed.  This  saves  both  storage  and  com¬ 
puting  time. 

Once  the  choice  of  geometries  to  be  investigated  has  been  made, 
the  flux  surfaces  are  generated  and  tessalated  in  a  semi-automatic 
fashion.  The  initial  physical  parameters  are  specified  and  the  system 
is  iterated  to  an  equilibrium  state  in  a  Lagrangian  fashion.  A  La- 
grangian  algorithm  was  chosen  for  its  simplicity  and  to  follow  the 
spatial  motion  of  the  flux  surfaces  under  different  initial  conditions. 

When  an  equilibrium  state  is  reached,  the  code  shifts  into  the 
Eulerian  portion  of  the  model  and  the  transport  and  diffusion  equations 
are  solved.  Diffusion  and  transport  involve  rapid  flow  along  the  mag¬ 
netic  surfaces,  especially  those  connected  to  the  limiter,  so  that  mass, 
energy  and  magnetic  flux  are  transported  across  the  cell  boundaries 
to  avoid  a  severe  distortion  of  the  grid  that  might  result  with  a 
Lagrangian  calculation. 

The  Lagrangian  force  balance  portion  of  the  model  has  been  dis¬ 
cussed  in  detail  in  Ref.  6,  Some  additional  general  remarks  are  made 
below.  This  report  is  primarily  concerned  with  the  Eulerian  transport 
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portion  of  this  model  and  the  application  of  this  model  to  the 
investigation  of  an  axisymmetric  poloidal  divertor  system. 

The  divertor's  function  is  to  remove  impurities  from  the  exterior 
of  a  magnetically  confined  plasma.  Divertor  efficiency  depends  on  the 
ability  of  the  plasma  in  the  scrape-off  layer  to  stream  into  the 
divertor  throat,  dumping  the  impurities  onto  the  limiter.  With  this 
model,  we  investigate  the  physics  of  the  scrape-off  layer  and  the 
limiter  by  calculating  the  plasma  flux  into  the  divertor  throat  and 
the  effect  of  charge-exchange  in  the  scrape-off  layer.  The  limiter 
must  be  able  to  remove  impurities, permanently.  Our  objective  is  to 
determine  the  factors  influencing  limiter  action,  and  the  importance 
of  features  such  as  a  neutral  blanket  on  the  operation  of  the  limiter. 

A  short  review  of  the  Lagrangian  portion  of  the  model  is  in  order. 
As  mentioned  above,  the  model  is  based  on  the  quasi-static  evolution 
of  force-balance  in  the  plasma.  The  Lagrangian  force-balance  is 
presently  calculated  with  only  the  pressure  gradient  and  the  j  x  B 
forces  (including  the  toroidal  and  poloidal  magnetic  fields) .  The 
mass,  the  entropy  and  the  magnetic  fluxes  are  all  conserved  as  the 
vertices  are  adjusted  in  position  in  response  to  any  force  imbalance. 

The  degree  of  adjustment  is  directly  proportional  to  the  degree  of  force 
imbalance.  This  adjustment  produces  a  negative  definite  change  in 
potential  energy  ensuring  the  approach  to  a  minimum  energy  state. 

The  initial  conditions  depend  on  the  problem  to  be  solved.  We 
usually  prescribe  the  initial  pressure,  density  and  temperature 
profiles.  These  functions  are  then  adiabatically  adjusted  as  the 


flux  surfaces  are  moved.  The  magnetic  fluxes  remain  conserved  during 
the  Lagrangian  motion. 

The  initial  pressure  distribution  is  prescribed  as  a  function  of 
the  flux  surface  variable  <)>.  The  pressure  must  remain  independent  of 
the  poloidal  variable  on  all  flux  surfaces  not  connected  to  tne 
limiter.  The  pressure  is  updated  adiabatically  (pV^  m  constant)  on 
these  flux  surfaces  where  is  the  volume  enclosed  by  the  flux 
surface.  The  pressure  is  updated  adiabatically  (pV^  =  constant)  on 
the  outer  surfaces  (those  connected  to  a  limiter)  vertex  by  vertex. 

Here  y_  is  the  volume  of  the  basic  cell  surrounding  the  vertex  i. 

The  toroidal  magnetic  field  is  initially  prescribed  as  a  function 
of  1/R  over  the  whole  system.  As  the  system  is  iterated  to  equilibrium, 
the  toroidal  flux  is  conserved  while  ensuring  RB0  =  constant  in  the 
vacuum  region  and  RBQ  =  F(ijj)  only  in  the  plasma.  The  ij>  dependence  of  F 
gives  the  deviation  of  the  toroidal  field  from  the  vacuum  field  and 
thus  measures  the  poloidal  current.  The  assumption  of  axisymmetry 
coupled  with  the  triangular  mesh  allow  the  currents  to  be  determined 
in  closed  form. 

The  initial  poloidal  flux  is  calculated  from  the  desired  toroidal 

current  distribution.  The  poloidal  flux  is  then  conserved  throughout 

the  iteration  to  force  balance.  Equilibrium  states  have  been  obtained 

for  both  low  and  high  poloidal  beta  discharges  and  for  doublet-like  and 

10 

and  poloidal  divertor  systems. 
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There  are  several  important  advantages  to  this  model  compared  to 
the  usual  MHD  equilibrium  model.  (1)  By  focussing  on  the  flux  sur¬ 
faces  we  avoid  the  difficulty  of  having  to  directly  solve  a  highly  non¬ 
linear,  elliptic  partial  differential  equation  with  complicated  boundary 
conditions.  (2)  Complicated  geometries  with  interfaces  and  separa- 
trices  can  be  accommodated  with  a  minimum  of  difficulty.  (3)  The 
poloidal  current  distribution  is  induced  in  a  self-consistent  manner 
so  as  to  ensure  force  balance.  (4)  The  Euler ian-Lagrangian  nature  of 
the  model  allows  both  the  gross  dynamics  of  the  plasma,  and  the  diffu¬ 
sion  and  transport  to  be  followed.  (5)  The  finite  resistivity  of  the 
plasma  and  shell  permit  the  simulation  of  discharges  for  times  exceeding 
the  penetration  time  of  magnetic  fields  through  the  shell. 

In  the  next  section  we  present  a  detailed  discussion  of  the  Euler- 
ian  transport  model.  These  equations  are  expressed  on  a  triangular 
mesh  in  Section  ill.  In  Section  IV,  we  discuss  time  scales  asso¬ 
ciated  with  the  problem  and  the  approximations  that  stem  from  the 
ordering  of  those  time  scales.  In  Section  V,  we  present  a  simple  analy¬ 
tical  model  involving  transport  along  the  flux  surfaces  connected  to 
the  limiter  (the  scrape-off  layer) .  Numerical  results  are  presented 
in  Section  VI,  with  a  summary  and  conclusions  in  the  last  section. 

II.  PHYSICAL  MODEL 
A.  Transport 

Four  fluid  equations  are  needed  to  describe  the  dynamics  of 
a  two-temperature  plasma.  The  continuity  equation  is 
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+  v  •  (pv)  =  o. 


(ii-i) 


The  momentum  equation  is 


P 


dv 

dt 


Vp+-jxB-F  -F  , 
e  c  —  —  n  — cx 


(II-2) 


where  F  is  the  viscous  drag,  and  F  is  the  ion  momentum  chanqe 
— r\  — cx 

resulting  from  charge  exchange.  The  energy  equations  are 
3T. 

n  — -  =  -  nV  •  (T.V.)  + ( 2-y)  nT.V  •  V.  -  V  •  (q. )  +  Q  -  Q_ 
3t  1-1 

3M  n 
e 


(II-3) 


(T 


Ti), 


and 


3T  3^ 

n  =  -  nV  •  (TV)  +(2-y)  nT  V  •  V  -  V  •  (q  )  +  — 
3t  e—e  e  ~e  «  ab 


( I X  —  4 ) 


.  3M  n 

+  —  -  j  •  V.  T  -  rr~  (T  -  T. )  , 


a 


b  e  M.T 
l  e 


where  O  is  the  heat  generated  by  the  viscous  drag,  Q  is  the  energy 

n  cx 

lost  by  charge  exchange,  the  symbol  b(i)  indicates  the  direction 
parallel  (perpendicular)  to  the  total  magnetic  field  and  is  the 

ion  (electron)  heat  flux. 
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Thu  viscous  dray ,  charge  exchange  and  heat  llux  terms  are  us  yet 
.ecj'i.ci  For  t  h**  calculations  which  follow,  wo  hov-  talon 


F  =  V  •  ir_  (II-S) 
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-  irv  T  -  K.1 7  T ,  -  Klb  X  7  T., 
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q  =  -  K'  V  T  -  K'  7  T  -  Y  'b  X  7  T  ; 

b  be.  ate  :  e 


(1'1-ii) 


r , 

■(,  k.4  f  K"  •  is  thu  classical  iori  (electron)  thermal  conductivity.  * 

V,  n  >  s  the  i  r.  it.  /  between  the  -ortmi  and  ionic  specie 

Ire  lotui,  ir.d  ee  owe  try  ar->  — >r>r>letely  determined  by  four  n.ote 


at  !.c'i'  s . 


One  is  ohm*;;  low 


including  the  e'ectron  pressure  and  tali 


{ IX— 9) 


E+-VxB=^  +  ^+  —  j  x  B  -  —  Vp  . 
—  c  —  —  o  o  nec  —  ne  e 

b  x 


We  need  three  of  Maxwell's  equations, 


V  •  B  =  0, 


(11-10) 


c  §  E  •  dl. 


(II-ll) 


§  B 


(H-12) 


where  we  use  the  integral  forms  of  Faraday's  and  Ampere's  laws.  These 
equations  are  solved  in  an  axisymmetric  orthogonal  flux  coordinate 
system  described  in  the  next  section. 

We  seek  to  solve  this  set  of  equations  on  the  resistive  diffusion 
timescale,  which  is  long  compared  to  the  Alfven  transit  time.  On  the 
resistive  time  scale  the  inertial  terms  in  the  momentum  equation  are 
negligible  and  the  system  evolves  through  a  series  of  quasi-static 
equilibrium  states  satisfying  the  static  force  balance  condition 

Vp  =  -  j  X  B.  (11-13) 

C  —  — 


The  method  of  attaining  an  equilibrium,  by  solving  Eq.  (11-13) ,  is 
discussed  in  detail  in  Ref.  6. 

Once  an  equilibrium  state  has  been  obtained,  the  transport  equations 
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are  solved  in  the  following  manner: 

(1)  The  poloidal  (toroidal)  current  is  obtained  by  integrating  Ampere's 
law  in  the  toroidal  (poloidal)  direction.  (2)  The  perpendicular  veloc¬ 
ity  is  obtained  from  Ohm's  law.  (3)  The  parallel  velocity  is  found  by 
assuming  steady-state  flow  along  the  flux  surface;  the  parallel  pressure 
gradient  is  balanced  by  parallel  viscous  drag  and  charge  exchange. 

(4)  The  electric  fields  are  determined  from  the  parallel  component  of 
Ohm's  law.  (5)  Once  the  currents,  velocities  and  electric  fields  are 
determined,  the  mass  density,  the  energies  and  the  magnetic  fluxes  are 
transported  according  to  Eqs.  II  -  1,  3,  4  and  11. 

The  preceding  algorithm  completely  determines  the  equilibrium  plasma 
properties.  The  algorithm  is  outlined  in  Figure  1.  In  the  next  section 
we  discuss  the  toroidal  coordinate  system  and  the  differential  algebra. 

B.  Coordinate  System 

We  employ  spherator-like 1 3 • 1 4  orthogonal  basis  vectors  determined 
by  the  poloidal  flux.  Axisymmetry  is  assumed  so  all  physical  quantities 
are  independent  of  the  toroidal  coordinate  (8).  Figure  2  illustrates 
these  basis  vectors.  Following  Ref,  14,  we  introduce  a  poloidal  flux 
function  f(R,Z,t]  defined  by 


(11-14) 


and  a  function  x(R/Z,t)  such  that 


Vf  •  Vx  =  o. 


(11-15) 
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That 


wit  a 


’.'hi  3 


is,  the  snri aces  of  constant  i  are  ortnogonal  to  the  surface 
ant  >. 

i:  magnetic  field  is  written  .a;. 


b 


p  +  b 
~A  ~~ 


Uv:  p^loiiaj  i  ioU  qjvon  by 


--  B,  +  r?i . 
x  i-  •- 


The  no  trie  of  the  f'i’jjv  coordinate  syten.  (  x»  f  >  ^ )  is  given  by 


(dl)?  =  vB  )'2(dv):‘  +  (RB  R?(d6);. 

>-  X 


gives  rise  to  the  following  differential  operations 


Vf  =  e,,  RB  -|  +  e  B  —■  +  §  l  ||  , 
V  X  <»V  X  X  "X  9  R  30 


2  3  2  3  1  3Ae 

-  Bx  W  B  ‘  +  BX  jx"  (B  '  +  k  30  ~  ' 

X  >■- 


Tl  >*r  BX  3 

7  *  A  -  *y  |r  3T  *  r  *  (kV 


1  A  T  !^0  i  9M 

J  *  §x  |_Bx  at  ~  R  39  J 


+  e  RB" 

-P 


A.,. 

—  ( — -)  -  -  (-A 

dX  RB  df  B 

X  X 


] 


dS  -  ey  |  dxdO  +  e  |  dYdfl  +  ev,  dxdt, 

X  X  X 


s  of 


(II-1B) 


(11-17) 


(11-18) 


(11-19) 


(11-20) 


(11-21) 


(11-22) 
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(11-23) 


dV  =  B-i  df'd  i.dt' 

A 


with  the  derivatives  given  by 


RB‘ 

X 


j 

if 


3  d 

B  —  -  B  — 
R  JZ  Z  JR 


(11-24) 


B2  — 

X  •>  X 


j _ 

3Z 


+  B 


R  3R 


(11-25) 


The  transport  equations  are  more  amenable  to  solution  in  a  coordi¬ 
nate  system  with  axes  parallel  and  perpendicular  to  the  magnetic  field. 
For  this  reason  we  introduce  a  coordinate  system  described  by 


and 


1 

h 


(e  +  fe  )  , 

0  A 


(11-26) 


e 


i' 


B 

where  f  =  ^  and  h2  =  1  +  f2.  e^  and  eg  define  the  plane  perpendicular 

to  the  magnetic  field  direction  . 

b 

C.  Currents 

We  divide  the  current  into  parallel  (j^)  and  perpendicular  (j^) 
parts  defined  by 


i 


13 


I 

.  A 


The  above  expressions  are  used  in  the  velocity  and  field  calculations. 

In  the  equilibrium  calculation  there  are  no  unresolved  forces  along  the 
tlux  surface,  so  the  f-component  of  the  current  density  is  zero. 

Using  the  mteqral  form  of  Amp4re's  law,  we  obtain  the  toroidal 
current  ( I  )  on  the  basic  cell  vertex.  All  vertex  functions  are  con- 

■i 

sidered  constant  over  the  basic  cell,  so  the  toroidal  current  density 
(}(  )  is  lust  the  current  divided  by  the  cell  area.  Applying  vnpere's 
law  to  the  toroidal  magnetic  field  and  integrating  around  the  torus,  we 
find  the  poloidal  current  flux  on  a  triangle  side.  The  K,Z  components 
of  the  current  density  are  found  on  the  triangle  in  a  manner  similar 
to  the  poloidal  magnetic  field  calculation.*5 

There  are  several  significant  features  of  these  algorithms.  The 
components  of  the  current  density  are  found  in  closed  form.  The  to¬ 
roidal  component  is  independent  of  the  toroidal  angle  and  thus  divergence 
free.  Summing  the  poloidal  current  fluxes  around  a  triangle  yields 
zero  assuring  that  no  current  sources  are  being  introduced,  and  thus 
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automatically  assuring  a  divergence-free  current.  The  electric  field 
calculation  is  discussed  in  the  next  section. 

D.  Electric  fields 

The  electric  fields  are  determined  from  the  parallel  component  of 
Ohm '  s  law 


Now 


1  - 
-  ev. 

ne  b 


(11-31) 


„  1  3A 

^  =  -  v*  -  c  aT 


(11-32) 


where  A  is  the  vector  potential  induced  by  the  ohmic  heating  transformer. 
Thus 


A  =  (0,0, Aq) . 


(11-33) 


The  magnetic  flux  is  contained  within  the  core  of  the  transformer, 
so  we  have 

B  =  V  X  A  =  0 

within  the  plasma  and  from  Eq.  (11-21) 

A.  =  A.  (R_1)  (11-34) 

W  u 
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Taking  the  induced  potential  to  be  a  linear  function  of  time,  we 


have  the  following  expression  for  the  toroidal  component  of  the  induced 
vector  potential 

A0  =  f  (11-35) 

where  a  is  a  constant  to  be  determined  from  the  ohmic  heating  current. 

The  toroidal  component  of  the  electric  field  (the  inductive  field) 
is  then 


E 


e 


(11-36) 


and  from  Eq.  (11-31)  we  have  the  expression  for  the  poloidal  electric 
field 


E 

X 


(VfV 


B  9p 
X  _ e  a 

ne  3X  fcR' 


(11-37) 


To  find  E^  we  need  an  expression  for  the  electrostatic  potential. 
Substituting  Eq.  (11-32)  into  Eq.  (11-31)  and  assuming  the  electron 
temperature  is  constant  on  a  flux  surface  (an  assumption  which  will  be 
relaxed  in  the  future)  we  find 


Letting 


V  (<t> 


T 


—  fnn) 
e 


$  = 


* 


T 

-  — —  inn 
e 


(11-38) 


(11-39) 


and  integrating  along  a  flux  surface,  we  obtain  the  electrostatic  poten¬ 
tial  to  within  a  surface  function. 15 


The  surface  function  ($m)  is  determined  by  requiring  the  absence  of 
any  net  plasma  rotation. 


ipv  ^  =  0. 

XBX 


(11-40) 


This  is  an  expedient,  which  allows  us  to  avoid  .excessive  complexity. 

The  physically  correct  determination  of  $  is  from  a  surface  integral 
over  the  component  of  the  force  equation.  12 

We  use  Eqs.  (11-46)  and  (11-47)  below  to  express  in  our  magnetic 
coordinate  system.  Then  applying  the  above  assumption,  we  have  the 
following  expression  for  the  derivative  of  the  surface  function 


where 


and 


=  -  (Dm  -  tr  I?  Amx')>/Gm,  (n-4D 


Rp  dx  _3 
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D(K)  = 


A(Y 


/  ch  B  I  b  hB0  I  nec  ne  x^eJJ 


■  «'>■/  » *<w. 


Gm  =  1-  £8— 


S(*)  =  f 


h2f  B  ' 


(11-42) 


(11-43) 


(11-44) 


The  ('-component  of  the  electric  field  is  given  by 


V^'X')  = 


RSX  *  (».X') 


RB^  ~  A(¥,X'>  -  RB ^  §j*¥). 


(11-45) 


E.  Velocities 

We  separate  the  velocities  into  a  parallel  and  perpendicular  part 
in  the  same  manner  as  the  currents,  thus 


v  =  7“  (V  +  fV  )  #  (11-46) 

b  h  6  X 

v  —  ~  (fV-  -  V  )  (11-47) 

sh6  x 

and 

V  =  V  +  V  (11-48) 

— x  —4'  -^s 

A  A  A 

With  these  choices,  (elfl,  e,  ,  e  )  is  a  right-handed,  orthonormal  basis. 

— b  — s 

We  shall  now  express  the  velocity  in  its  components  on  this  basis. 

We  obtain  the  perpendicular  components  from  Ohm’s  law. 


V 
~ 1 


XB  + 


—  B  X  2  B2 


-  B  X  Vp 

ne  -  el 


(11-49) 


This  expression  has  two  components  which  are 


and 


B0h2 

\  ’  S^h 


(11-50) 


(11-51) 


Thus  is  determined  when  E  and  p  are  known. 

The  parallel  velocity  is  obtained  from  the  parallel  component  of 
the  momentum  equation.  Dotting  Eq.  (II-2)  with  efa  and  assuming  the 
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acceleration  term  is  negligible,  i.e.,  the  parallel  velocity  has 
attained  its  terminal  value,  we  have 


6b  *  VP  = 


F  +  e 
— n  b 


F 

~cx 


(11-52) 


or  that  the  parallel  pressure  gradient  is  balanced  by  the  viscous  drag 
and  charge  exchange  terms. 

The  full  viscous  stress  tensor  for  this  coordinate  system  is  pre¬ 
sented  in  the  Appendix.  The  parallel  component  of  the  viscous  drag  term 
is 


(F  K  =  K 

n  b  b 


f  _2  3  /  ”bb\ 

nb  h  BX  9X  (  B  )  ' 


(II-53) 


where  it.  is  the  e,  £  component  of  the  stress  tensor  and  n,  is  the  paral- 
bb  b  b  b 

lei  viscosity.  Keeping  the  dominant  term  in  the  stress  tensor,  we  obtain 


(F  =  4  ^  4  b2  T" 

nb  b  3  b  h  x  3  X 


(11-54) 


The  charge  exchange  term  describes  the  momentum  change  due  to  a  charge 
exchange  event.  It  can  be  expressed  as  the  product  of  a  coefficient  and 
the  relative  velocity  between  the  ion  and  neutral  species.11  The 
parallel  component  of  this  momentum  change  is 


(F  ).  =  n  mo  VjVl 
cx  b  o  cx—b1  — 'b 


(11-55) 
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where  m  is  the  reduced  mass  and  o  is  a  function  of  the  relative  kinetic 

cx 

11  12 

energy.  > 

Equation  (11-52)  then  takes  the  form 


b 


n  mo  V  I  V  I 
o  cx  — b  — b 


(11-56) 


which  is  solved  on  each  flux  surface  for  the  parallel  velocity  on  a 
triangle. 

With  the  velocities  determined  on  the  triangles,  the  mass  density 
is  advanced  in  time  on  the  vertices,  through  use  of  the  continuity  equa¬ 
tion.  The  ion  and  electron  energies  are  time-advanced  using  Eqs.  (II-3) 
and  ( 1 1—4 ) .  The  integral  form  of  Faraday's  law  (Eq.  (11-11))  is  used  to 
find  the  new  magnetic  flux.  At  this  point  all  quantities  have  been 
advanced  to  the  new  time.  A  flow  chart  for  the  complete  algorithm  is 
illustrated  in  Figure  3.  The  computational  model  for  these  equations 
is  developed  in  Section  III. 

III.  COMPUTATIONAL  MODEL 

A.  Introduction 

The  set  of  equations  developed  in  Section  II  will  be  approximated 
by  finite  differences  on  a  triangular  mesh.  The  variables  in  these 
equations  will  be  represented  as  triangle,  side,  or  vertex  quantities 
on  this  mesh.  This  differencing  procedure  is  somewhat  complex.  We 
will  illustrate  it  by  discussing  some  important  basic  concepts.  We 
begin  with  the  concept  of  the  basic  computational  cell.  Next  we  define 
the  gradient  operator,  and  finally  we  will  demonstrate  the  finite 
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differencing  of  the  nonlinear  diffusion  equation.  Further  information 
on  triangular  gridding  can  be  found  in  Refs.  6,  16,  17. 

The  basic  computational  cell  is  the  shaded  region  shown  in  Fig.  4. 
It  is  formed  by  joining  the  side  bisectors  of  the  triangles  surrounding 
the  general  vertex.  Any  physical  quantity  defined  on  the  vertex  is 
considered  to  be  constant  throughout  the  basic  cell,  and  each  triangle 
surrounding  the  general  vertex  contributes  1/3  of  its  area  to  the  area 
of  the  basic  cell. 

We  now  illustrate  how  a  gradient  is  represented  in  this  model. 

If  vertex  quantities  are  linear  functions  of  position,  then,  given  the 

function  g  (defined  on  vertex  m) ,  the  function  g  at  any  other  point, 
m 

n,  say,  can  be  written,  without  approximation,  as 


9n  =  9m  +  R  *  v9n 
n  m  ti  n 


(III-l) 


Here  R  is  the  vector  from  the  location  of  g  to  the  chosen  point. 

— n  m 

Now  consider  the  triangle  j  defined  by  two  side  vectors  S^, 
with  the  vertex-defined  quantities  g,  g^  and  g^+^.  The  index  j  indi¬ 
cates  triangle  quantities,  the  index  i  vertex  quantities.  (See  Fig.  5.) 
Following  Ref.  16,  the  gradient  of  g,  uniquely  defined  on  the  triangle 
and  constant  throughout,  is 


Vgj  =  (grg)  h+i  -  (gi+rg>  5i 


(HI-2) 


=  i  9< 


n  X  (r' 

_ -1-1  —1+1 

2A  . 
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+ 

where  £  is  the  side  vector  S  rotated  clockwise  by  n/2  radians,  and  is 
shorthand  for  a  cross  product  with  a  basis  vector.  Here  n  is  a  unit 
vector  normal  to  the  computational  plane  and  is  the  area  of  the 
triangle.  The  conclusion  is  that  gradients  may  be  naturally  represen¬ 
ted  on  triangles,  and  easily  calculated  on  them. 

Now  consider  the  nonlinear  diffusion  equation 


i£  _ 

at 


v 


(X  Vg)  +  H, 


(HI-3) 


where  g  is  a  vertex  quantity  and  X  and  H  (the  source  function)  are 
assumed  to  be  constant  throughout  a  triangle.  Again  following  Ref.  16, 
we  introduce  the  flux  of  the  diffusing  quantity  within  each  triangle  as 


F\  =  -  X_.Vg_..  (III-4) 

Gauss's  theorem  says  the  volume  integral  of  the  left  side  of  Eq.  ( X 1 1— 3 ) 
over  the  basic  cell  is  equal  to  the  integral  of  the  normal  component  of 
F  over  the  boundary  of  the  basic  cell.  The  flux  contribution  G_.  from 
triangle  j  (Fig.  6)  is  given  by 


G  . 
3 


<S+  -S+. 

—  l+l  —  : 


)  ■ 


(III-5) 


Summing  around  the  central  vertex,  the  finite  difference  form  of 
Eq.  ( III-3)  is 
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At  n  Ls=l  S  1  S  j  =  l  ^  3 

y  a  -*  — ^ 


(in-6) 


3=1  3 


where  the  coupling  coefficient 


“S  =  hlX3+l  Ctn9j+1  +  Aj-1  ctnej-i^ 


and  the  angles  0..^  lie  opposite  the  side  S  and  a_.  =  1/3  A^ .  Examina¬ 


tion  of  Fia.  6  will  reveal  that  all  these  quantities  are  area  or  side 


subelements  of  a  basic  cell.  Expressions  for  the  divergence,  curl  and 


Lap lac i an  operators  are  presented  in  Ref.  17.  These  are  all  the 


difference  operators  we  need.  We  now  proceed  to  represent  the  equa¬ 


tions  of  Section  II  on  a  triangular  grid. 


B.  Currents 


With  the  magnetic  fields  given  on  the  mesh,  the  currents  can  be 


calculated  directly  from  them  by  using  Ampere's  law. 


The  poloidal  magnetic  field  is  assumed  constant  throughout  each 


triangle.  Integrating  it  around  the  perimeter  of  a  basic  cell  gives 


the  toroidal  current  (I  K  at  the  corresponding  vertex.  As  a  basic 


cell  (vertex)  function,  IQ  is  constant  over  the  cell.  Thus  the 

U 


toroidal  current  density  (defined  on  a  triangle)  is  given  by 


3  1  Ie/Ai' 

1=1  1 


(III-7) 
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The  poloidal  current  is  found  on  triangle  sides,  by  integrating 
the  toroidal  magnetic  field  around  the  torus.  The  Lagrangian  portion 
of  the  code  assures  the  physical  condition  that 

RB  =  £(?) . 

d 

Thus  the  dependence  of  RB„  is  a  measure  of  any  induced  poloidal  cur- 

p 

rents.  It  gives  the  deviation  of  the  toroidal  field  from  the  vacuum 
field.  To  make  this  quantitative,  apply  Ampdre's  law  to  a  triangle 
side  connecting  two  flux  surfaces.  The  poloidal  current  flux  on  that 
side  is 


I  =  ~  [R.i./A.  - 

y  4n  111 
s 


i+1 


°i+l/Ai+lJ 


(III-8) 


where  is  the  toroidal  magnetic  flux  at  vertex  i. 

I  is  thus  the  poloidal  current  flux  through  a  ribbon  obtained 
*s 

by  revolving  a  line  of  length  ds  between  the  V  and  V  +  dS1  flux 
surfaces  about  the  axis  of  symmetry.  This  is  equivalent  to  the  manner 
in  which  the  poloidal  magnetic  flux  is  calculated5  and  thus  the  poloidal 
current  density  can  be  determined  on  the  triangles  from  the  relation 
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2nr  e9 


X  VI  . 
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(III-9) 


Again  the  current  density  is  an  average  of  currents,  divided  by  a 
triangle  area. 

The  assumption  of  axisymmetry  allows  the  current  densities  to  be 
calculated  in  closed  form  on  the  triangular  grid.  The  toroidal  current 
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density  is  independent  of  the  toroidal  coordinate  and  automatically 
divergence-free.  Since  the  poloidal  current  is  the  same  for  all  sides 
interconnecting  two  flux  surfaces,  summing  the  currents  around  the 
triangle  yields  zero,  assuring  that  no  current  sources  are  being  intro¬ 
duced.  There  are  no  unresolved  forces  along  the  flux  surface  in  the 
Lagrangian  calculation,  so  the  ¥  component  of  the  current  density  is 
zero. 

Expressions  for  the  electric  field  are  presented  in  the  next 
section. 


C.  Electric  Fields 

Equations  (11-36)  and  (11-37)  yield  the  toroidal  and  poloidal 
electric  fields  as  triangle  quantities 

Ee  =  -  £  (III-10) 

and 

1  1  /  ^Pe  ^Pe\  E0 

Ex  =  (je+fV  -  7ZT  [br  IT  +  Bz  if)  '  T  (III'U) 

The  procedure  for  finding  E^  is  somewhat  more  involved.  All  the 
integral  quantities  (Eq.  11-41)  are  weighted  to  the  vertices  giving 
the  electrostatic  potential  on  the  vertex.  The  gradient  of  the  poten¬ 
tial  then  yields  E  on  the  triangles. 

D.  Velocities 

The  perpendicular  components  of  the  velocities  are  calculated  as 
triangle  quantities  using  Eqs.  (11-50)  and  (11-51) , 
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The  parallel  component  of  the  velocity  is  found  by  assuming  steady- 
state  flow  along  the  flux  surface.  The  parallel  pressure  gradient  is 
then  balanced  by  the  parallel  viscous  drag  and  charge  exchange  terms. 
This  calculation  is  discussed  in  greater  detail  in  Sections  IV  and  V. 

To  facilitate  the  gradient  and  divergence  operations,  the  ion 
(fluid)  and  electron  velocities  are  decomposed  into  R-  and  Z-components 
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We  now  have  all  the  tools  we  need  to  solve  the  transport  and 
diffusion  equations. 


E.  Mass  and  Energy  Transport 

As  noted  previously,  the  transport  and  magnetic  diffusion  are 
Eulerian.  The  triangles  are  held  fixed  in  position  and  the  mass. 
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energies  and  magnetic  fluxes  are  transported  among  the  cells. 

The  continuity  equation  (Eq.  II-l)  is  differenced  in  the  same 
manner  as  the  diffusion  equation  of  Section  IIIA.  The  new  vertex 
mass  density  is 


n  n 

pn  =  P°  +  —  At  A'  p  .  (vAr'-v  Az'  )  ./  A  a 
i  1  2  .  ,  i  Z  R 

3=1 


3  1=1-  3 


( IX  I- 14 ) 


where  (AR7  AZ7  )  are  the  components  of  the  side  of  triangle  j  opposite 
the  vertex  i. 

The  energy  equations  (Eqs.  II-3,  II-U)  are  solved  in  a  similar 
fashion.  We  include  both  the  parallel  and  perpendicular  viscous  terms. 
The  perpendicular  term  will  be  important  (because  of  the  long  path 
length)  in  the  scrape-off  layer.18 
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Finite-differencing  Eq.  (II-3)  ,  we  have  new  vertex  ion  temperature 
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where 


(III-17) 


(III-18) 


(III-19) 


and  K1  and  n .  can  be  found  in  Section  II. 

1 

The  magnetic  flux  is  integrated  forward  in  time  using  the  integral 
form  of  AmpSre's  law 


:^E 


dfc 


( III-21) 


where 


$  =  /  B  •  dA 


The  new  toroidal  electric  field  is  defined  on  the  vertices.  Inte¬ 
grating  this  field  around  the  torus,  the  new  poloidal  magnetic  flux 
through  a  triangle  side  is 
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(II 1—22 ) 


tpn  =  ^J°  +  2lt  cAt  ( R .  E 


i“0  .  -  Vl  E0.  > 

1  1+1 


where  i,  i+1  are  the  vertices  at  the  ends  of  side  s.  The  V  •  B  =  0 

condition  is  assured  since  the  sum  of  the  fluxes  around  a  triangle 
5,6 

is  zero. 

We  wish  to  keep  the  toroidal  flux  a  vertex  quantity,  so  Eq. 
(III-21)  is  integrated  around  the  basic  cell  to  give  the  new  toroidal 
flux 


.n 

$  .  =  <t> .  + 

l  l 


Ate 
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(III-23) 


Now  that  the  mass,  energy  and  magnetic  flux  have  been  trans¬ 
ported,  the  system  is  no  longer  in  equilibrium  so  we  return  the 
Lagrangian  portion  of  the  code  to  iterate  to  a  new  equilibrium  state. 
The  new  transport  coefficients  are  found,  and  the  code  reenters  the 
transport  portion  of  the  model. 


F.  Exactness  of  Conservation  Equations 

The  basic  computational  cell  is  irregular  (see  Fig.  7)  .  Therefore 
one  might  question  whether  conservation  equations  can  be  accurately  re¬ 
presented  by  finite  difference  operations  on  these  cells.  This  section 
shows  that  accurate  representations  are  used  in  the  present  model. 

Probably  the  simplest  numerical  question  is  the  order  of  accuracy 
of  a  difference  algorithm  on  a  triangular  grid.  Fritts1^  has  demon¬ 
strated  second-order  accuracy  for  a  wave  equation.  We  use  his  methods 
in  this  code. 

A  second  natural  question  is  whether  "conserved  quantities"  are 
conserved  to  machine  roundoff.  If  this  property  is  desired  (and  the 
authors  believe  it  is  in  this  case)  ,  it  can  easily  be  accomplished  by 
expressing  the  changes  in  the  conserved  quantities  in  terms  of  fluxes 
between  the  cells,  and  insuring  that  the  flux  is  calculated  in  the  same 
way  from  both  "sides"  at  a  cell  boundary.  This  is  done  in  the  code. 

There  is  another  important  aspect  to  "conservation"  in  a  magnetic 
geometry.  It  is  conservation  within  a  flux  surface.  In  a  tokamak-like 
system,  the  velocity  components  in  the  magnetic  surfaces  may  be  two  or 
three  orders  of  magnitude  larger  than  the  component  perpendicular  to 
the  surface.  Thus  a  numerical  algorithm  may  seriously  degrade  conser¬ 
vation  of  energy  (or  density)  in  a  flux  surface,  if  it  allows  even  a 
small  part  of  the  surface  velocity  to  appear  as  a  normal  velocity  in 
the  different  equations.  The  remainder  of  this  section  explains  how 
this  problem  is  avoided  in  a  triangular  gridding  system. 

Consider  a  Gedanken  experiment  where  the  temperature  is  zero 


everywhere  except  at  one  vertex.  There  it  is  T'.  A  set  of  flux  sur¬ 


faces  is  illustrated  in  Figure  7a.  We  will  show  that  our  Eulerian 
scheme,  for  advancing  the  energy  and  mass  densities  forward  in  time, 
preserves  the  physics  of  the  conservation  equations;  i.e.,  that  the 
parallel  and  perpendicular  transports  are  completely  uncoupled. 

In  our  Gedanken  experiment  we  shall  consider  parallel  heat  flux  in 
the  energy  equation,  and  show  that  no  energy  is  transported  between  flux 
surfaces.  The  energy  equation  has  the  form 


3T 

at 


=  v 


=  -  v 


(^b 
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(III-23) 


where  is  the  parallel  conductivity  tensor.  From  Eq.  (III-2  ),  the 
temperature  gradient  is  zero  everywhere  except  for  the  triangles  around 
vertex  "a"  where  it  is  normal  to  the  sides  opposite  vertex  "a". 

Finite  differencing  Eq.  III-23  on  vertex  "c",  we  have 
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( III-24) 


For  triangles  1.  and  4.,  the  parallel  heat  flux  is  zero  since  the  temp¬ 
erature  gradient  is  normal  to  the  flux  surface.  For  triangle  2.,  the 
parallel  heat  flux  is  normal  to  the  surface  area  vector  of  the  basic 
cell.  Thus  if  we  are  solving  for  the  parallel  transport,  we  are  assured 
that  no  energy  is  being  transported  across  the  flux  surface  numerically 
The  mass  continuity  equation  also  exactly  conserves  density.  Here 


we  assume  the  density  is  zero  everywhere  except  at  one  vertex  and  assume 
only  flow  in  the  magnetic  surface.  See  Figure  7b.  The  continuity 


equation  has  the  form 


_L_ 
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E  P  .  v* (nxS  ) 
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( III-25) 


where  (the  density  on  the  triangle)  is  the  average  of  the  densities 
at  the  vertices  of  the  triangle.  From  Eq.  11-56,  we  see  that  the  paral¬ 
lel  velocity  is  a  function  of  the  parallel  component  of  the  pressure 
gradient.  For  triangles  1.  and  4.,  the  parallel  velocity  is  zero.  The 
parallel  flux  of  triangle  2  is  again  normal  to  the  surface  of  basic 
cell  c,  and  no  mass  is  tranported  across  the  flux  surface. 

Finally,  we  may  assure  that  if  we  have  only  perpendicular  flow,  no 
mass  is  transported  along  a  flux  surface.  Considering  only  perpendicu¬ 
lar  flow,  the  change  in  the  mass  density  at  vertex  b  is 


^b 

At 


where  p. 

3 

Now 


1/3  p'.  See  Figure  8. 


v 


v  e  +  v  e 
f  y  s  s 


and  the  velocity  arising  directly  from  the  pressure  is 

*  M  ®x  *  VP 

and 

vs  =  hue^  •  Vp 

(c.f.  Eqs .  11-50,  51)  where 


( I 11—26 ) 
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the  perpendicular  velocity  can  then  be  expressed  as 


v  =  P  S  X  Vp. 
— x  0 


(III-28) 


Letting  F .  be  the  flux  contribution  from  each  of  the  triangles 
surrounding  vertex  b,  we  have  for  triangles  2  and  3  respectively 

F2  =  3  P'  M  ?P2  ‘  h  (HI-29) 

and 

F3  =  J  p'  M  Vp3  •  ( III-30) 

with  the  pressure  gradients  given  by 


and 


Vp2  =  1T2  e6  x  *3 


Vp3  =  iJT  60  *  ^5 


(III-31) 


(III-32) 


A  and  A  are  the  triangle  areas, 
2  3 


and 


*2  ’  2  «e  •  ‘^3  x  h' 


A3  =  2  ®9  •  (54  *  — S* ' 


Substituting  into  Equation  III-26  and  summing  around  vertex  b 
we  have 
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which  is  identically  zero. 


r 


A  similar  argument  applies  to  the  perpendicular  heat  flux  since 
the  conductivity  tensor  contains  both  a  ip  and  an  s  component. 

The  power  of  the  triangular  gridding  scheme  is  immediately  evident. 
Not  only  can  complicated  geometries  be  modified  with  a  high  degree  of 
resolution,  but  the  gradient  operator,  as  a  triangle  function,  is  exact 
and  there  is  absolutely  no  numerical  coupling  between  the  parallel  and 
perpendicular  transport. 

G.  Summary 

We  now  summarize  the  transport  calculation.  The  toroidal  and  poloi- 
dal  components  of  the  current  density  are  found  as  triangle  quantities 
from  Ampere's  law. 

The  poloidal  and  toroidal  electric  fields  are  found  as  triangle 
quantities  from  the  parallel  component  of  Ohm’s  law. 

Once  the  parallel  fields  and  current  densities  are  known,  the  perpen¬ 
dicular  components  of  the  velocity  are  determined  as  triangle  quantities 
from  the  perpendicular  component  of  Ohm's  law.  The  parallel  component  is 
determined  as  a  triangle  quantity  by  equating  the  parallel  viscous  drag 
and  charge  exchange  with  the  parallel  pressure  gradient. 

The  perpendicular  electric  field  is  then  found  by  solving  the  paral¬ 
lel  component  of  Ohm's  law  for  the  electrostatic  potential  under  the  con¬ 
dition  that  there  is  no  net  plasma  rotation. 

Now  that  the  dynamical  variables  have  been  determined, the  transport 
equations  are  solved  to  give  the  new  mass  densities  and  energies  on  the 
vertices  and  the  new  toroidal  (poloidal)  flux  as  a  vertex  (side)  quantity. 

In  the  next  section  we  discuss  a  reduced  set  of  equations  that  stem 
from  an  examination  of  the  time  scales  associated  with  this  problem. 
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IV.  SIMPLIFIED  MODEL 


A.  Time  Scales 

Our  model  is  based  on  the  assumption  that  the  equations  pre¬ 
sented  in  the  last  two  sections  can  be  solved  on  the  resistive  time 

scale  (t  ).  Since  the  resistive  time  scale  is  much  larger  that  the 

o 

Alfven  (t^)  time  scale,  the  inertial  terms  in  the  equation  of  motion 
are  negligible  and  the  system  can  be  taken  as  evolving  through  a  se- 
quenceof  equilibrium  states  satisfying  the  force-balance  condition.  By 
comparing  the  characteristic  times  for  the  other  dissipative  processes 
with  the  Alfven  transit  time,  we  find  that  density  and  temperature  per¬ 
turbations  on  a  magnetic  flux  surface  cannot  be  maintained  for  times  on 
the  order  of  the  Alfv6n  transit  time. 

The  Alfven  transit  time  is  given  by 


T 

A 


(B2/4irp}1/2 


(IV-1) 


where  a  is  the  plasma  radius.  The  characteristic  times  for  the  impor¬ 
tant  dissipative  processes  are: 

Resistive  diffusion. 


t  ~  l2/d 

a  B  2 

Particle  diffusion, 
L2m. 

i 

T  rs*  . 

D  t.T.  ' 

1  l 


(IV-2) 


(IV- 3) 
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Thermal  conduction. 


T 

K 


L^n 


(IV-4) 


and  for  viscous  drag 


T 

n 


(iv-d) 


where  L  is  a  characteristic  scale  length  and  i  is  the  ion-ion  colli- 

i 

sion  time. 

Comparing  the  resistive  diffusion  time  to  the  Alfven  transit  time, 
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where  g  =  nT/  (B2/8tt)  ,  S2.  =  eB/m.c,p  =  v  /Q.  and  v  is  the  ion-ion 

i  l  e  the  e  i 

collision  frequency.  The  criterion  required  for  our  basic  assumption 
to  be  met  is  thus  well  satisfied. 

Comparing  the  characteristic  times  of  the  other  dissipative  pro¬ 
cesses  along  the  flux  surfaces  to  the  Alfven  time,  we  have 


td/ta 


ta/ta  ~ 


Tn/TA  ~  ^  aT  <<  x- 
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The  above  implies  that  density,  temperature  and  velocity  perturba¬ 
tions  cannot  exist  on  a  flux  surface  for  times  on  the  order  of  the 
Alfven  transit  time;  thus,  to  a  very  good  approximation  the  density  and 
temperature  can  be  taken  to  be  constant  on  a  flux  surface.  This  applies 


to  all  flux  surfaces  except  for  the  separatrix  and  the  scrape-off  layer 
where  the  parallel  transport  must  be  determined.  This  problem  is 
addressed  in  Sections  IV. C  and  V. 

B.  Perpendicular  Transport 

As  an  examination  of  the  time  scales  indicates, the  particle  and 
heat  flux  is  very  rapid  along  the  flux  surfaces  compared  to  the  resis¬ 
tive  diffusion  across  the  surfaces.  To  lowest  order  then,  we  need  only 
solve  for  the  perpendicular  transport  and  take  the  parallel  transport  as 
occurring  with  infinite  speed. 

Our  simplified  model  then  consists  of  including  only  the  perpendic¬ 
ular  transport  coefficients  (cf .  Eq.  II-8) .  Once  the  continuity  and 
energy  equations  have  been  solved,  the  density  and  temperatures  are 
averaged  on  the  flux  surfaces.  The  only  modifications  to  the  equations 
presented  in  the  previous  two  sections  are  that  all  the  parallel  trans¬ 
port  coefficients  are  set  to  zero. 

In  the  next  section  we  discuss  the  parallel  transport  in  the  scrape- 
off  layer  where  this  approximation  cannot  be  made. 

C.  Simplified  Transport  Equations 

In  the  interior  of  the  discharge,  the  ordering  discussed  in  the 
previous  section  insure  that  temperature  and  pressure  are  constant  on 
each  magnetic  surface.  As  one  proceeds  outward  toward  the  limiter,  the 
effect  of  the  external  boundaries  is  to  increase  transport  in  the  sur¬ 
faces,  and  perturb  this  equilibrium. 

The  pressure  and  temperature  distributions  on  the  outer  surfaces 
are  still  nearly  constant,  but  they  exhibit  some  variation  on  these 
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surfaces,  accompanied  by  very  rapid  transport  on  the  (resistive)  times 


of  interest  to  us.  Thus  it  is  appropriate  to  calculate  the  transport 
in  the  outer  surfaces  by  determining  the  fluxes  across  magnetic  sur¬ 
faces  from  the  averaged  quantities,  and  to  determine  the  transport  in 
the  surfaces  from  the  quasi-stationary  solutions  of  the  equations  in¬ 
cluding  surface  variations. 

We  shall  assume  that  this  outer  region  is  cold  enough  that 
viscosity  is  no  longer  important,  and  that  the  neutral  density  is  suf¬ 
ficiently  high  that  Ohmic  heating  is  negligible.  Then  equations  (II-l 
to  4)  reduce  to 


V  •  (nv)  =  S, 


(IV-8) 


Vp  -  —  j_  x  B  +  nrnv/icx  =  0, 


( IV-9) 
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Here  t  is  the  charge-exchange  time  and  S  represents  a  source  term, 
cx 

Ohm's  law  plays  an  important  role  in  this  region,  since  an  electro¬ 
static  potential  is  required  near  the  limiter  to  enforce  ambipolarity . 

To  be  consistent  with  the  preceding  assumptions,  we  must  neglect  the 
conductivity  and  ohmic  heating  vector  potential  terms.  We  also  neglect 
the  electron  pressure  and  Hall  terms.  Then  Eq.  (II-9)  becomes 
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with  <f>  the  electrostatic  potential. 

With  B  determined  separately  and  using  the  zeroth  order  current 
densities  as  determined  from  Ampdre's  law,  equations (IV-8) through (IV-12 ) 
plus  charge  conservation 

V  •  j_  =  0  ( IV-13 ) 

can  be  solved  for  the  first  order  quantities  p,  T^,  T.,  ,  v  and  j_. 

This  solution  is  discussed  in  the  next  section. 

V.  AN  ANALYTIC  BENCHMARK 

A.  Ordering 

The  solutions  of  the  transport  equations  in  the  exterior  will  in¬ 
clude  rapid  flow  and  heat  flux  along  B  toward  the  limiter.  Little  else 
is  easy  to  discern  from  the  transport  equations.  Thus  it  is  appropriate 
to  reduce  these  equations  to  an  analytically  tractable  form.  This  can 
most  easily  be  done  by  assuming  only  small  variations  on  the  surfaces. 

We  assume  that  n  and  T  satisfy  inequalities  like 
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where  y  is  a  measure  of  length  in  the  poloidal  direction  and  #dy  =  L 

For  simplicity  we  shall  also  assume  T  =  T.  =  T,  and  examine  only  one 

e  i 

temperature  equation.  Then  the  equations  we  wish  to  solve  are 
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where  w  =  e^  x  v  x  e^  represents  the  two-velocity  components  in  the 
magnetic  surface.  We  now  examine  them  order  by  order. 

Represent  each  of  these  nine  physical  variables  by  a  series  in  €, 
for  example, 

n  =  n(0>  +  £n(1)  •••  (V-7) 

Assume  the  source  terms  are  first  order  in  (E.  Further,  assume  that  the 
zero-order  density  and  temperature  are  surface  functions.  Then  equations 
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plus  an  energy  equation  which  vanishes  identically. 

These  equations  may  be  solved  as  follows.  Equation  (V-8)  may  be 
integrated,  yielding  a  velocity  vector  with  a  prescribed  poloidal  varia¬ 
tion,  plus  a  constant  vector  in  the  ignorable  direction.  The  parallel 
component  of  Eq.  (V-10)  requires  the  parallel  component  of  w  to  vanish. 
The  only  consistent  solution  is 

wl0)  =  0.  (V-12) 

By  similar  analysis,  the  integral  of  Eq.  (V-9)  and  the  perpendicular 
component  of  Eq.  (V-10)  uniquely  determine 
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and  finally,  from  Eq.  (11), 


(0) 
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This  completely  determines  the  lowest  order  solution. 

B.  First-Order  Equations 

In  this  ordering,  the  flow  and  surface  variations  of  density  and 
temperature  appear  in  first  order.  We  shall  represent  the  density  and 
temperature  source  terms  through  sources  S  and  q  which  are  given. 
Then  the  first-order  equations  are 
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We  take  and  j  ^  as  given  quantities  as  determined  from 

0  X 

Ampere's  law.  The  previous  set  of  equations  is  then  solved  in  the 
following  manner. 

From  Eq. (V-22) 
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The  electrostatic  potential  is  found  from  Eq.  (V-21) , 
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The  other  velocity  components  are  then  found  from  Eqs.  (V-15)  and  (V-20) , 
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The  H*  and  x  components  of  the  first  order  current  densities  and  given 
by  Eqs.  (V-19)  and  (V-10) , 
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The  momentum  equation  then  yields  the  pressure 
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and  the  toroidal  current  density 
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The  remaining  equation  to  solve  is  the  energy  equation.  Writing  p 
as  (n  nT  +  n  nT  ) ,  it  can  now  be  expressed  in  terms  of  known 
quantities . 

In  other  calculations  of  this  nature,  consistency  conditions 
determine  the  integration  constants  in  the  velocity  and  current.  Here 
they  are  determined  by  the  sheath  conditions  at  the  limiter. 

C.  Boundary  Conditions 

Boundary  conditions  are  required  for  equations  (V-24) ,  (V-25) , 
(V-28)  and  (V-29) .  These  boundary  conditions  are  determined  by  the 
limiter . 

Integrating  Eq.  (V-24) , 


<P(X) 


dX  +  4>£ 


(V-31) 


where  is  the  sheath  potential  at  the  limiter.  This  potential  will 
generally  be  a  few  Debye  lengths  (Aq)  thick  and  is  of  the  form 


~  4irne 


and  thus  only  a  function  of  the  temperature  at  the  limiter. 


In  a  similar  manner,  Eq.  (V-25)  gives 
/ 

X  /  — 

V  (/>  «  B  (/)  f  ( 

x  *  J  xA 


(0). 


*X0)b2x  n 


(0)  3Y 


/nRV""\  B  (x'> 

l-Tn*’- b— 


(V-32) 


where  v^ 


is  the  ion-thermal  velocity  at  the  limiter. 
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VI.  NUMERICAL  RESULTS 
A.  Equilibrium 

The  initial  flux  surfaces  for  the  divertor  grid  are  generated 
by  a  multicurrent  vacuum  code  with  the  positions  and  strengths  of 

external  magnetic  field  coils  chosen  to  be  representative  of  PDX  II 

19 

geometry.  Seventeen  surfaces  are  used  in  the  grid  (each  current 
axis,  defined  by  a  single  point,  is  counted  as  a  separate  surface). 
Each  surface  is  then  divided  into  a  number  of  straight-line  segments, 
the  end  points  of  which  are  the  triangle  vertex  positions.  The  com¬ 
puter  code  then  automatically  generates  the  triangular  grid  structure 
from  these  vertex  positions.  The  initial  grid  is  shown  in  Figure  9. 

The  two  outermost  surfaces  represent  the  conducting  shell 
(copper) .  Proceeding  inward,  the  next  surface  is  a  vacuum  surface, 
the  next  two  surfaces  comprise  the  scrape-off  layer,  the  next  surface 
is  the  separatrix  and  the  inner  surfaces  describe  either  the  plasma 
discharge  or  the  divertor. 

The  power  of  the  triangular  gridding  scheme  in  providing  a 
mesh  with  variable  resolution  is  immediately  evident.  Flow  in  the 
scrape-off  layer,  the  region  bounded  by  the  separatrix  and  the  plasma- 
vacuum  interface,  is  expected  to  be  quite  rapid  and  thus  this  region 
should  be  finely  resolved.  This  region  is  divided  into  382  triangles. 
Flow  in  the  main  body  of  the  plasma  discharge,  the  area  enclosed  by 
the  separatrix  and  to  the  right  of  the  null  point,  should  be  much 
less  than  the  flow  in  the  scrape-off  layer  and  thus  the  region  does 


not  need  to  be  so  finely  resolved.  The  plasma  discharge  region, 
which  occupies  a  much  larger  area  than  the  scrape-off  layer  is 


represented  by  373  triangles.  In  all,  the  divertor  grid  is  comjjosed 


of  17  flux  surfaces,  777  vertices,  1458  triangles  and  2234  line 
segments . 

Once  the  grid  structure  has  been  generated,  the  plasma 
profiles  are  calculated  with  the  pressure,  density , temperatures  and 
plasma  current  specified  on  the  vertices.  The  initial  current  density 
and  temperature  profiles  are  parabolic  out  to  the  separatrix.  The 
density  profile  has  a  cubic  radial  dependence.  These  profiles  then 
fall  off  nearly  exponentially  from  the  separatrix  to  small  but  non¬ 
zero  values  at  the  plasma-vacuum  interface.  The  initial  density  and 
temperatures  are: 

n(r  =0)  =  9  x  10^  cm  ^ 


and 


T  (r  =  0)  =  T.  (r  =  0)  =  0.7  KeV. 
e  l 

At  the  plasma-vacuum  interface 

,  ,  12  -3 

n  =  1.7  x  10  cm 

and 


T  =  T.  =  0.01  KeV. 
e  l 

The  system  is  then  iterated  to  an  equilibrium  state  as 
discussed  in  Section  I  of  this  report  and  in  detail  in  Ref.  6.  The 
toroidal  magnetic  field  is  slowly  ramped  over  10  iteration  steps  to 
its  final  value  of  20  KG.  Depending  on  the  initial  plasma  parameters, 
the  system  attains  equilibrium  after  100-150  iteration  steps. 
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The  equilibrium  grid  structure  is  presented  in  Figure  10  with 

the  corresponding  profiles  in  Fiqures  11  —  1  1 .  The  plasma  current  has  a 

value  of  990  KA  and  the  central  density  and  temperatures  are 

n(r  =  0)  =  9.3  x  10^  cm  ^ 

and  T  (r  =  0)  =  T. (r  =  0)  =  .8  KeV. 

e  i 

The  major  radius  is  123  cm  and  the  minor  radius  is  60  cm.  The 
poloidal  beta  is  0.7.  The  width  of  the  scrape-off  layer  (at  the  null 
point)  on  the  outside  of  the  separatrix  is  16  cm  and  the  width  on 
the  inside  is  22  cm.  The  perimeter  of  the  plasma  discharge  has  a 
length  of  370  cm  (see  Figure  14).  The  results  of  the  transport 
calculation  are  presented  in  the  next  section. 
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B.  Transport 


Now  that  an  equilibrium  state  has  been  attained,  the  diffusion 
and  transport  equations  can  be  solved.  Before  presenting  the  results 
of  this  calculation  several  remarks  are  in  order. 

As  noted  in  Section  I,  one-dimensional  transport  calculations 
which  necessitate  the  determination  of  flux  surface  averaged  quantities 
are  not  appropriate  for  investigating  the  divertor  system.  This  is 
because  the  surface  inteqrals  generally  involve  terms  of  the  form 
1  and  the  poloidal  field  (B^)  goes  to  zero  at  the  null  point  on 
the  separatrix.  Without  a  careful  choice  of  computational  grid  and 
metric,  most  surface  averaged  quantities  are  formally  infinite  on  the 
separatrix . 

This  difficulty  is  avoided  completely  in  our  computation  for 

several  reasons.  First,  the  calculation  is  fully  two-dimensional  and 

surface  averaged  quantities  are  not  required.  Second,  the  poloidal 

magnetic  field  is  calculated  as  a  triangle  function  (constant  over  the 

area  of  the  triangle),  and  although  B  is  zero  at  the  x-point  it  is 

X 

always  non-zero  on  the  trianqles.  As  a  result,  all  the  flow  velocities 

(also  triangle  functions)  are  finite.  Third,  from  the  metric  used 

in  this  calculation  (Eq.  11-18)  the  quantity  dy/B  is  just  the  lenqth 

X 

(in  centimeters)  along  a  flux  surface;  and  thus,  the  numerical 
integration  required  to  obtain  the  electrostatic  potential  (Eq.  11-41) 
is  finite.  Simply  put,  the  computational  algorithm  treats  the  x-point 
in  the  same  way  as  any  other  point  in  the  system. 


4  H 


s 


Purely  classical  transport  coefficients  are  used  in  the 
calculations  presented  here.  Ideally,  neoclassical  coefficients 
should  be  used,  especially  in  the  scrape-off  layer  but  this  is  not 
expected  to  drastically  affect  the  nature  of  the  results.  The 
boundary  conditions  on  the  limiter  are  chosen  to  be  the  simplest 
possible.  The  density,  temperature  and  poloidal  velocity  go  to 
zero.  Neoclassical  coefficients  and  more  realistic  boundary  condi¬ 
tions  on  the  limiter  will  be  incorporated  in  any  future  modifications 
of  the  code. 

One  of  the  most  important  questions  concerning  a  discharge 
that  has  a  separatrix  is  the  role  of  the  separatrix  in  determining 
the  diffusive  flow  and,  in  particular,  does  the  separatrix  present 
any  barrier  to  the  flow.  Our  results  indicate  that  the  separatrix 
does  not  hamper  the  diffusive  flow  across  the  flux  surface.  This 
result  is  obvious  by  examing  the  expression  for  the  flow  across  a 
flux  surface  (Eq.  11-50).  Rewriting  Equation  11-50,  we  have 


and  it  is  clear  that  is  not  zero  on  the  separatrix,  contrary  to 

results  obtained  from  surface-averaged  calculations. 

This  conclusion  is  borne  out  in  the  numerical  results  as 

shown  in  Figure  15.  Here  we  have  plotted  the  vectors  for 

approximately  half  the  triangles  just  interior  to  the  separatrix. 

2  2 

The  maximum  r (z)  -  component  of  is  6.12  x  10  cm/sec  (6.16  x  10 
cm/ sec) .  This  is  the  length  of  the  "unit"  vector  and  all  vector 
lengths  are  scaled  to  this  "unit"  length.  The  diffusive  velocity 


across  the  separatrix  is  nearly  uniform,  around  the  major  portion 

of  the  separatrix  except  near  the  null  point  where  it  is  somewhat 

larger  than  the  average  since  f  (f  =  B  /B. )  is  small  there.  V„ 

X  6  t 

takes  on  values  between  200  and  700  cm/sec  with  the  larger  values 
near  the  null  point. 

The  x -component  of  the  diffusion  velocity  in  the  scrape- 

off  layer  is  shown  in  Figure  16.  The  magnitude  of  V  ranges 

X 

3  5 

between  5  x  10  cm/sec  and  1  x  10  cm/sec  with  the  larger  values 
appearing  near  the  null  point.  Note  that  the  stagnation  point  is 
far  removed  from  the  limiter.  This  gives  rise  to  a  very  long  flow 
path  on  the  outside  of  the  plasma  discharge.  Also  note  that  the 
particle  flow  rate  is  somewhat  larger  on  the  outside  of  the  discharge. 
This  increased  flow  rate  is  probably  due  to  the  smaller  toroidal 
field  on  the  outside.  This  result  may  be  somewhat  fortuitous 
physically  in  that  the  long  flow  path  on  the  outside  may  be  compen¬ 
sated  for  by  the  larger  velocity  there.  This  phenomenon  is  being 
investigated  further  by  incorporating  a  more  realistic  neutral  density 
profile  in  the  scrape-off  layer  and  more  realistic  boundary  conditions 
on  the  limiter. 

A  measure  of  divertor  efficiency  can  be  made  by  comparing 

the  diffusive  flow  across  the  separatrix  to  the  diffusive  flow  into 

the  divertor  throat.  The  average  velocity  across  the  separatrix  is 

2 

~  4  x  10  cm/sec. 

The  length  of  the  separatrix  is 

2 

L  =3.70x10  cm. 
s 
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Thus  a  measure  of  the  diffusive  flux  across  the  separatrix  can  be 
written  as 


D  =  <VU)>  L 
s  r  s 

~  5  2 

~  1.5  x  10  cm  /sec. 

The  average  flow  velocity  along  the  scrape-off  layer  is 

/V  «  5  x  104  cm/sec, 

X 

and  taking  the  width  of  the  divertor  throat  to  be 


W  ss  16  cm, 
d 

the  average  flux  into  the  divertor  throat  is 


D 


d 


<V  >  W 
X 


5  2 

~  8  x  10  cm  /sec. 

The  the  divertor  is  capable  of  accepting  the  plasma  deposited  in  the 
scrape-off  layer. 

The  pressure,  density  and  temperature  profiles  after  several 
diffusion  times  are  illustrated  in  Figures  17-19. 
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VII.  SUMMARY  AND  CONCLUSIONS 

Ke  have  developed  a  fully  two-dimensional  equilibrium  and 
transport  computer  simulation  model  and  applied  it  to  the  investi¬ 
gation  of  the  poloidal  divertor  system.  The  model  is  unique  in  that 
it  uses  as  a  finite  difference  mesh  a  general  connectivity  triangular 
grid.  The  advantages  of  this  model  are  numerous. 

1)  Complicated  geometries  can  be  handled  with  a  minimum  of 
difficulty. 

2)  The  resolution  across  the  mesh  is  highly  variable. 

3)  The  Eulerian-Lagrangian  nature  of  the  model  allows  the 
gross  dynamics  of  the  plasma  and  diffusion  and  transport  to  be 
followed. 

4)  The  finite  resistivity  of  the  plasma  and  shell  permit 
the  simulation  of  discharges  for  times  exceeding  the  field 
penetration  time  of  the  shell. 

5)  The  poloidal  currents  are  induced  in  a  self-consistent 
manner  so  as  to  ensure  force  balance. 

6)  The  currents  are  found  in  closed  form. 

7)  The  gradient  operator,  as  a  triangle  function,  is  exact. 

8)  In  spite  of  the  orders  of  magnitude  difference  between 
the  flow  along  and  across  the  flux  surfaces,  the  conservation 
equations  exactly  conserve  mass  and  energy. 

9)  Numerically,  there  is  nothing  peculiar  about  the  x-point. 

Applying  this  model  to  an  investigation  of  the  poloidal 

divertor  system  we  are  able  to  show  that  physically  there  is  also 
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nothing  peculiar  about  the  x-point;  that  is,  the  major  portion  of 
the  plasma  discharge  is  not  affected  by  the  presence  of  an  x-point. 
The  diffusive  flow  across  the  separatrix  is  nearly  uniform,  around 
the  perimeter  of  the  separatrix.  We  have  found  that  there  is  a 
very  long  flow  path  on  the  outside  of  the  discharge.  This  may  be 
compensated  for  somewhat  by  the  increased  poloidal  velocity  in  this 
region.  The  flux  across  the  separatrix  is  comparable  to  the  flux 
into  the  divertor  throat.  A  conclusion  as  to  whether  the  divertor 
can  remove  impurities  effectively  and  permanently  must  await 
modification  to  the  code  to  include  various  neutral  density 
profiles  and  a  sheath  potential  on  the  limiter. 


' 

! 


ACKNOWLEDGMENTS 


Drs.  John  L.  Johnson  and  Michio  Okabayashi  of  the 
Princeton  Plasma  Physics  Laboratory  have  contributed  numerous  valuable 
suggestions  to  the  development  of  this  model.  Ufa  Christiansen 
provided  magnetic  field  calculation  routines  which  speeded  its 
development.  Conversations  with  Drs.  Martin  Fritts  and  Jay  Boris 
at  NRL  were  especially  useful.  It  is  a  pleasure  to  acknowledge 
other  useful  discussions  with  Dale  Mead,  John  Greene,  Ray  Grimm, 

Alan  Boozer,  Steve  Jardin,  Joan  Ogden,  Teruyo  Tamano,  and  Abe  Kadish. 
This  work  was  supported  by  the  U.S.  Department  of  Energy  under 
Contract  EX-76-A-36-1006. 


54 


APPENDIX 


In  this  section  we  present  the  components  of  the  viscous  stress 
tensor  in  our  coordinate  system.  In  the  momentum  equation  we  need  only 
the  parallel  component  of  the  viscous  drag  but  we  must  include  both  the 
parallel  and  perpendicular  components  in  the  ion  energy  equation.  18 
The  parallel  viscous  drag  term  has  the  form21 


F 

— b 


4/3 


n  e  v2v 

b  b  b  b 


(A-l) 


which  can  be  expressed  as 


4  =  4/3  Vb  h  Bx  k  (f  ST* 


(A-2) 


where  is  the  parallel  viscosity  coefficient, 8 


n,  =  0.96  n .  T .  t  . 

b  ill 


and  i.  is  the  ion  collision  time. 

l 

In  the  e  ,  e  ,  e  basis,  the  components  of  the  viscous  stress 
b  s  T 

tensor  are2  1 
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where 


the  velocity  is 


,  1  niTi 

n  =  I  — ' 

i 


v  =  v  e  +  v,  e,  +  v  e 
—  V  f  b  b  s  s 


(A— 4) 


and 
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Expanding  equations  A-3,  the  components  of  the  viscous  stress  tensor  are 


Tn  =  -  %s  =  n'[\  hLJiZnl£‘)  "  VsRBX^''  ‘  <Bx>  ^ 
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The  total  heat  generated  by  tne  viscosity  is  given  by 


Q  =  if  :  Vv. 
n  ~ 


(A-7 ) 
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Expanding  the  above,  we  have 
rfB  3v  \  2 


Q  =  V3  n 


(fB  3v  \ 
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This  expression  is  quite  complicated  and  an  ordering  must  be 
imposed  to  make  it  more  tractable.  Physically  we  know 

Vv  <<:  V  V 

And,  in  general,  the  logarithmic  derivatives  are  smaller  than  the  other 
terms  in  the  expression  since  they  involve  "geometric"  factors. 

We  calculate  the  total  viscous  stress  tensor  according  to 
Equation  (A-3)  but  subsequently  drop  all  the  logarithmic  terms  in 
Equation  (A-5)  to  find  an  approximate  expression  for  the  viscous  heating 
term 


Q  »  4/3  n, 
n  b 


B  9v  ' 
h  9x 


( A-9 ) 


This  completes  the  terms  we  need  to  include  viscosity  in  our  model. 
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Fig.  2  -  The  toroidal  coordinate  system,  e*  is  along  the  flux  surface  in 
the  direction  of  the  poloidal  magnetic  field,  e  is  perpendicular  to  the 
flux  surface,  e  is  in  the  toroidal  (ignorable^direction.  The  poloidal 
field  is  measured  from  the  outside,  and  (e  ,  e  ,  e  )  form  a  right-handed 
coordinate  system.  *  b 


Define  Profiles 


Fig.  3  -  Flow  chart  of  complete  computational  algorithm.  The  code  is  split 
into  a  Lagrangian  part  (force-balance)  and  an  Eulerian  part  (transport). 
The  Lagrangian  iteration  takes  place  within  the  complete  timestep  loop. 
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Fig.  4  -  Detail  of  triangular-grid  elements.  A  triangle  is  made  up  of 
three  directed  line  segments,  and  a  vertex  represents  the  (shaded)  region 
defined  by  the  center  of  mass  of  each  surrounding  triangle  and  the  midpoint 
of  each  side. 
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Fig.  6  -  Solving  the  diffusion  equation  on  a  triangular  mesh.  The  basic 
cell  function  (g  )  is  integrated  forward  in  time  by  integrating  the  triangle 
functions  (XVg)  Ground  the  perimeter  of  the  basic  cell,  m^  and  m  are  the 
portions  of  the  side  bisectors  of  triangle  j  which  contribute  to  the  boundary 
of  the  basic  cell. 


Fig.  7a  -  Parallel  thermal  flux  calculation.  Temperature  is  assumed 
non-zero  only  at  vertex  a.  The  arrows  represent  AT'.  The  parallel 
thermal  fluxes  of  triangles  1  and  4  are  zero.  The  parallel  thermal 
flux  of  triangle  2  is  normal  to  the  surface  of  basic  cell  c.  As  a 
result,  no  thermal  energy  is  transported  across  the  flux  surface. 


Fig.  7b  -  Parallel  mass  flux  calculation.  The  mass  density  is  assumed 
non-zero  only  at  vertex  a.  The  parallel  velocity  is  directly  proportional 
to  the  parallel  pressure  gradient.  The  arrows  represent  vj,.  The  parallel 
mass  fluxes  of  triangles  1  and  4  are  zero.  The  parallel  mass  flux  of 
triangle  2  is  normal  to  the  surface  of  basic  cell  c.  As  a  result  no  mass 
density  is  transported  across  the  flux  surface. 
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Fig.  8  -  Perpendicular  mass  flux  calculation.  The  mass  density  is  assumed 
non-zero  only  at  vertex  a.  For  vertex  b,  the  perpendicular  flux  contribution 
from  triangle  2  is  cancelled  exactly  by  the  perpendicular  flux  contribution 
from  triangle  3.  As  a  result,  no  mass  density  is  transported  along  the  flux 
surface.  Note  the  counterclockwise  direction  of  the  side  vectors  which  define 
the  triangles. 
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Fig.  9  -  Initial  divertor  grid  configuration.  The  ?  direction,  pointing 
away  from  the  axis  of  symmetry,  is  "up".  The  two  outermost  surfaces  make 
up  the  copper  shell.  Moving  inward:  The  next  surface  is  a  vacuum  surface, 
the  next  twn  surfaces  comprise  the  scrape-off  layer,  the  next  surface 
is  the  separatrix.  The  other  surfaces  make  up  the  divertor,  on  the  left, 
and  the  plasma  discharge,  on  the  right.  Note  the  high  resolution  in  the 
scrape-off  layer. 
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Ftp,.  10  -  Equilibrium  grid  configuration.  The  outer  surfaces  have  moved 
away  from  the  wall  giving  a  clearer  picture  of  the  grid  structure.  Note 
that  the  null  point  is  well  defined. 
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Fig.  11  -  Pressure  profile.  A  two-dimensional  perspective  plot  of  the  pressure 
The  point  of  view  is  from  the  outside  of  the  discharge;  i.e.,  large  r,  looking 
into  the  divertor.  The  r,  z  dimensions  are  in  centimeters. 
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Fig.  16  -  Parallel  flow.  The  vectors  represent  the  x-component  of 
the  diffusive  velocity  in  the  scrape-off  layer.  The  vector  lengthy 
are  scaled  to  the  unit  length  which  has  components  (VR  =  9.71  x  10 
cm/sec,  V  =  1.04  x  10^  cm/sec).  The  stagnation  region  is  far  from 
the  limiter  giving  rise  to  a  very  long  flow  path  on  the  outside  of 
discharge . 


Fig.  17  -  Pressure  profile.  A  two-dimensional  perspective  plot  of  the 
pressure  after  several  diffusion  times.  Note  the  somewhat  flatter  profile 
across  the  main  body  of  the  discharge. 
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Fig.  20  -  Pressure  profile.  A  two-dimensional  perspective  plot  of  the 
pressure  as  viewed  from  the  inside  of  the  discharge,  i.e.  small  r. 
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